// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

#include <avtDistanceToBestFitLineExpression.h>
#include <vtkDataArray.h>
#include <avtDataTree.h>
#include <avtParallel.h>
#include <math.h>
#include <avtExecutionManager.h>

#define N_SUM   0
#define X_SUM   1
#define Y_SUM   2
#define XY_SUM  3
#define X2_SUM  4

#define N_CALC_VALUES 5

// ****************************************************************************
// Method: avtDistanceToBestFitLineExpression::avtDistanceToBestFitLineExpression
//
// Purpose: 
//   Constructor for the avtDistanceToBestFitLineExpression class.
//
// Arguments:
//   v : If true then difference vertically - otherwise calculate perpendicular
//       distance.
//
// Programmer: Brad Whitlock
// Creation:   Fri Nov 18 16:15:57 PST 2005
//
// Modifications:
//   
// ****************************************************************************

avtDistanceToBestFitLineExpression::avtDistanceToBestFitLineExpression(bool v) :
    avtBinaryMathExpression()
{
    verticalDifference = v;
    pass = 1;
    sums[N_SUM] = 0.;
    sums[X_SUM] = 0.;
    sums[Y_SUM] = 0.;
    sums[XY_SUM] = 0.;
    sums[X2_SUM] = 0.;
}

// ****************************************************************************
// Method: avtDistanceToBestFitLineExpression::~avtDistanceToBestFitLineExpression
//
// Purpose: 
//   Destructor for the avtDistanceToBestFitLineExpression class.
//
// Programmer: Brad Whitlock
// Creation:   Fri Nov 18 16:16:49 PST 2005
//
// Modifications:
//   
// ****************************************************************************

avtDistanceToBestFitLineExpression::~avtDistanceToBestFitLineExpression()
{
}

// ****************************************************************************
// Method: avtDistanceToBestFitLineExpression::PreExecute
//
// Purpose: 
//   Sets all of the sums to zero.
//
// Programmer: Brad Whitlock
// Creation:   Fri Nov 18 16:17:07 PST 2005
//
// Modifications:
//    Jeremy Meredith, Thu Feb 15 11:55:03 EST 2007
//    Call inherited PreExecute before everything else.
//   
// ****************************************************************************

void
avtDistanceToBestFitLineExpression::PreExecute(void)
{
    avtBinaryMathExpression::PreExecute();

    sums[N_SUM] = 0.;
    sums[X_SUM] = 0.;
    sums[Y_SUM] = 0.;
    sums[XY_SUM] = 0.;
    sums[X2_SUM] = 0.;
}

// ****************************************************************************
// Method: avtDistanceToBestFitLineExpression::Execute
//
// Purpose: 
//   Performs a multi-pass execute. In pass 1 we calculate sums for the best
//   fit line. In pass 2, we calculate the difference between the variables
//   and the best fit line.
//
// Programmer: Brad Whitlock
// Creation:   Fri Nov 18 16:17:24 PST 2005
//
// Modifications:
//   
//  David Camp, Thu May 23 12:52:53 PDT 2013
//   Changed the avtSIMODataTreeIterator::Execute call to support the new
//   threading changes. Now the data tree is not returned, but passed in as an
//   argument. Also you need to make sure all execution has completed before
//   continuing. This occured twice in this function.
// ****************************************************************************

void
avtDistanceToBestFitLineExpression::Execute(void)
{
    //
    // Sum the values required to do the best fit line.
    //
    pass = 1;
    avtDataTree_p tree    = GetInputDataTree();
    totalNodes = 2 * tree->GetNumberOfLeaves();
    avtDataTree_p newTree;
    avtSIMODataTreeIterator::Execute(tree, newTree);
    avtSIMODataTreeIterator::FinishExecute();
    newTree = 0;

    // Sum the array values over all processors, making sure each processor
    // gets the results
    double d[N_CALC_VALUES];
    SumDoubleArrayAcrossAllProcessors(sums, d, N_CALC_VALUES);
    for(int i = 0; i < N_CALC_VALUES; ++i)
        sums[i] = d[i];

    //
    // Make it perform the expression.
    //
    pass = 2;
    avtDataTree_p newTree2;
    avtSIMODataTreeIterator::Execute(tree, newTree2);
    avtSIMODataTreeIterator::FinishExecute();

    SetOutputDataTree(newTree2);
}

// ****************************************************************************
// Method: avtDistanceToBestFitLineExpression::DoOperation
//
// Purpose: 
//   Does the work of calculating the expression.
//
// Arguments:
//   in1 : The data array containing the X coordinate.
//   in2 : The data array containing the Y coordinate.
//   out : The resulting data.
//   ncomps : The number of components.
//   ntuples : The number of tuples.
//
// Note:       Note that we only populate the out data array in pass 2. The 
//             results of pass 1 get thrown away.
//
// Programmer: Brad Whitlock
// Creation:   Fri Nov 18 16:18:09 PST 2005
//
// Modifications:
//   Kathleen Biagas, Wed Apr 4 12:10:10 PDT 2012
//   Change float to double.
//
//   Alister Maguire, Tue Nov  8 12:44:55 PST 2016
//   Added thread_sums for thread safety and added
//   mutex locks where appropriate. 
//
// ****************************************************************************

void
avtDistanceToBestFitLineExpression::DoOperation(vtkDataArray *in1, 
    vtkDataArray *in2, vtkDataArray *out, int ncomps, int ntuples)
{

    // thread_sums is used to calculate the sum values
    // before adding them to sums. This is needed when
    // threading is enabled. 
    double thread_sums[N_CALC_VALUES]; 
    for (int i = 0; i < N_CALC_VALUES; ++i)
        thread_sums[i] = 0.;
   
    if(pass == 1)
    {
        // Sum up the values required to calculate the best fit line.
        thread_sums[N_SUM] += double(ntuples);
        for(vtkIdType i = 0; i < ntuples; ++i)
        {
            double x = in1->GetTuple1(i);
            double y = in2->GetTuple1(i);

            thread_sums[X_SUM] += double(x);
            thread_sums[Y_SUM] += double(y);
            thread_sums[XY_SUM] += double(x * y);
            thread_sums[X2_SUM] += double(x * x);
        }
    
        VisitMutexLock("avtDistanceToBestFitLineExpression::DoOperation");
        for (int i = 0; i < N_CALC_VALUES; ++i)
             sums[i] += thread_sums[i];
        VisitMutexUnlock("avtDistanceToBestFitLineExpression::DoOperation");
    }
    else if(pass == 2)
    {
        double dY = (sums[N_SUM] * sums[XY_SUM] - sums[X_SUM] * sums[Y_SUM]);
        double dX = (sums[N_SUM] * sums[X2_SUM] - sums[X_SUM] * sums[X_SUM]);
        double m, b;

        if(dX == 0.)
        {
            double xLine = sums[X_SUM] / sums[N_SUM];
            for(vtkIdType i = 0; i < ntuples; ++i)
            {
                double x = in1->GetTuple1(i);
                out->SetTuple1(i, x - xLine);
            }
        }
        else
        {
            m =  dY / dX;
            b = (sums[Y_SUM] - m * sums[X_SUM]) / sums[N_SUM];

            if(verticalDifference)
            {
                for(vtkIdType i = 0; i < ntuples; ++i)
                {
                    double x = in1->GetTuple1(i);
                    double y = in2->GetTuple1(i);
                    double yLine = m * x + b;
                    out->SetTuple1(i, y - yLine);
                }
            }
            else // perpendicular distance.
            {
                for(vtkIdType i = 0; i < ntuples; ++i)
                {
                    double x0 = in1->GetTuple1(i);
                    double y0 = in2->GetTuple1(i);
                    
                    // We know line eq. Y = MX + B
                    // Reorganized: ax + by + c = 0 where a=M, b=-1, c=B
                    //
                    // Dist from point to line:
                    //   |a*x0 + b*y0 + c| / sqrt(a^2 + b^2)
                    //
                    // Subst a,b,c
                    //   |M*x0 -1*y0 + B| / sqrt(M^2 + 1)

                    double d = fabs(m * x0 - y0 + b) / sqrt(m*m + 1);
                    out->SetTuple1(i, d);
                }
            }
        }
    }
}
